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Abstract 

Theoretical approaches to strong held phenomena driven by plasmonic helds 
are based on the length gauge formulation of the laser-matter coupling. From 
the theoretical viewpoint it is known there exists no preferable gauge and 
consequently the predictions and outcomes should be independent of this 
choice. The use of the length gauge is mainly due to the fact that the quantity 
obtained from hnite elements simulations of plasmonic helds is the plasmonic 
enhanced laser electric field rather than the laser vector potential. In this 
paper we develop, from hrst principles, the velocity gauge formulation of 
the problem and we apply it to the high-order harmonic generation (HHG) 
in atoms. A comparison to the results obtained with the length gauge is 
made. It is analytically and numerically demonstrated that both gauges 
give equivalent descriptions of the emitted HHG spectra resulting from the 
interaction of a spatially inhomogeneous held and the single active electron 
(SAE) model of the helium atom. We discuss, however, advantages and 
disadvantages of using diherent gauges in terms of numerical efficiency. 

Keywords: Strong held phenomena, time dependent Schrodinger equation, 
plasmonic helds 


1. Introduction 

Nowadays there exists a high demand for coherent light sources extending 
from the ultraviolet (UV) to the extreme ultraviolet (XUV) spectral ranges. 
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These sources provide important tools for basic research, material science 
and biology among other branches [1]. An important obstacle preventing 
these sources from reaching high efficiency and large duty cycles is their 
demanding infrastructure. The recent demonstration of XUV generation 
driven by surface plasmon resonances, conceived as light enhancers, could 
provide a plausible solution to this problem |2]. The high-order-harmonic 
generation (HHG) in atoms using plasmonics fields, generated starting from 
tailored metal nanostructures, requires no extra amplihcation of the incoming 
pulse. By exploiting the so-called surface plasmon polaritons (SPP), the 
local electric fields can be enhanced by several orders of magnitude Mi, 
thus exceeding the threshold laser intensity for HHG generation in noble 
gases. One additional advantage is that the pulse repetition rate remains 
unaltered without any extra pumping or cavity attachment. Furthermore, 
the high-harmonics radiation, generated from each nanostructure typically in 
the UV to XUV range, acts as a source with point-like properties, enabling 
collimation or focusing of this coherent radiation by means of constructive 
interference. This opens a wide range of possibilities to spatially arrange 
nanostructures to enhance or shape the spectral and spatial properties of the 
outgoing coherent radiation in numerous ways. 

One can shortly describe the high-order-harmonic generation based on 
plasmonics fields as follows (a more exhaustive description can be found in 
the seminal paper of Kim et ah 0): a femtosecond low-intensity laser pulse 
is coupled to the plasmon mode of a metal nanostructure inducing a col¬ 
lective oscillation of the free electrons within the metal. These free charges 
redistribute the electric field of the laser around each of the nanostructures, 
thereby forming a spot of highly enhanced electric held, also known as hot 
spot. The plasmon amplihed held exceeds the threshold of HHG, thus by 
injection of a gas jet, typically a noble gas, onto the spot of the enhanced 
held, high order harmonics from the gas atoms are generated. In the original 
experiment of Kim et ah [2], the output laser beam emitted from a low- 
power femtosecond oscillator was directly focused onto a 10 x 10 pm size 
array of bow-tie nanoantennas with a pulse intensity of the order of 10^^ 
W/cm^, which is about two orders of magnitude smaller than the threshold 
intensity to generate HHG in noble gas atoms. The experimental result of 
Ref. [2] showed that the held intensity enhancement factor exceeded 20 dB, 
i.e. the enhanced laser intensity is two orders of magnitude larger than the 
input one, which is enough to produce from the 7th to the 21st harmon¬ 
ics of the fundamental frequency by injecting xenon gas. For the case of 
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the laser wavelength corresponding to a Ti:Sa laser, i.e. about 800 nm, the 
wavelength of the emitted coherent radiation is between 38 nm and 114 nm. 
Additionally, each bow-tie nanostructure acts as a point-like source, thus a 
three-dimensional (3D) arrangement of bow-ties should enable us to perform 
control of the properties of generated harmonics, e.g. their polarization, in 
various ways by exploiting interference effects. Due to the strong conhnement 
of the plasmonic hot spots, which are of nanometer size, the laser electric held 
is clearly no longer spatially homogeneous in this tiny region. Since typically 
electron excursions are of the same order as the size of this region, important 
changes in the laser-matter processes occur, see. e.g. EIEIE]. 

So far, all of the the numerical approaches to study laser-matter processes 
in atoms and molecules driven by plasmonic helds, in particular HHG and 
ATI, are based on the length gauge of the laser-coupling formulation El El 
[IOl[IIl[l2l[l3l[Il[l5l[l6l[I71[l8l[l9l[20l|2llE2l[2a The use of the length 
gauge is mainly due to the fact that the quantity obtained from hnite elements 
simulations of plasmonic helds is the plasmonic enhanced laser electric field 
rather than the laser vector potential. Only a couple of papers employed an 
extension of the Strong Field Approximation (SEA), where an approximate 
version of the velocity gauge was used Diherent descriptions of light 

matter interaction (c.f. Ref. [28]), which include the full spatial dependence 
of the electromagnetic held, are closely related to the problem presented 
in our contribution. There are, however, distinct diherences amongst the 
general formulation of the non-dipole problem with the one we will tackle in 
the present article. For instance, the next order of the non-dipole description 
includes both the electric quadrople and the magnetic dipole terms, which 
are not present in our plasmonic helds, because the typical laser intensities 
are far below the ones needed to consider relevant these ehects. 

In this article, we concentrate our ehort on the formulation and numerical 
implementation of the velocity gauge description of light-matter interaction 
driven by plasmonic helds. From a pure theoretical viewpoint, it is known 
the velocity gauge is more appropriate and consequently our contribution 
will hll the missing gap, completing the whole picture in the modeling of 
laser-matter processes driven by plasmonic helds. 

The paper is organized as follows. In Sec. II, we shall present the velocity 
gauge formulation of the problem and we relate it to the length gauge, clearly 
showing the compatibility between them. The numerical implementation is 
presented in Sec. Ill, joint with a set of examples and a discussion about how 
the two diherent algorithms, i.e. the spectra split operator and the Crank- 
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Nicolson, behave as a function of the relevant parameters. Furthermore an 
analysis of the computational efficiency and scaling of both formulations is 
presented here. The paper ends with a short summary and an outlook. 

2. Theory and gauge transformation 

Quantum mechanics governs the evolution of the systems, atoms and 
molecules in our case, when they interact with an extental electromagnetic 
held. In particular, the Time Dependent Schrodinger Equation (TDSE) [33] 
allows us to obtain the complete time-space evolution of the particles. From 
a mathematical viewpoint, there are two different, but equivalent, expres¬ 
sions for the Hamiltonian which describes the interactions of the whole sys¬ 
tem. As a consequence the laser-matter problem can be formulated both in 
the so called velocity gauge (VG) or in the length gauge (LG), indistinctly. 
Formally, both gauges present equivalent descriptions of the quantum prob¬ 
lem [33], and therefore the results should not change if either the VG or LG 
is utilized to compute the observables of interest. Here, we detail how the 
gauge transformation is commonly implemented in the laser-matter interac¬ 
tion and in particularly when a spatial inhomogeneous held interacts with 
an atomic or molecular target. In general, we are interested in to describe 
the electron dynamics of an atomic or molecular system when it interacts 
with an electromagnetic held. For this case the TDSE reads (atomic units 
are used throughout the paper otherwise stated): 

HT(r,t) = z^T(r,t), (1) 

where, H, is the Hamiltonian of the quantum system and T(r,t) is the elec¬ 
tron wavefunction (EWF). 

Let us dehne the Hamiltonian, Hy, in the minimum coupling or VG for 
the electromagnetic held-matter interaction as: 

Hy = ^ [p + A(r,f)]VHo(r), (2) 

where, p = —iV, denotes the canonical momentum operator, A(r,t), is the 
vector potential of the electromagnetic held, which in this case corresponds 
to a spatial inhomogeneous or plasmonic held. In Eq. ([^, Vo(r,t) is the 
electrostatic Goulomb interaction between the charged particles. The vector 
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potential for the spatial inhomogeneous held typically can be represented in 
the following form: 


A(r,f) = [1+ e^(r)]Ah(f), 

Ah(f) = Aof{t) sm{uot + (pcEp)ez- 


( 3 ) 


Here, Ah(f), denotes the homogeneous or conventional vector potential, Aq 
is the amplitude of the vector potential, coq, is the central frequency, (pcEP 
is the carrier-envelope-phase (CEP) and f{t) is a function which dehnes the 
time envelope of the held, e is a small parameter that governs the strength of 
the spatial inhomogeneity (see e.g. [5] for more details) and g{r) describe the 
spatial dependence of the plasmonic held. Note that in the limit when e = 0, 
the vector potential held does not depend on the spatial coordinate anymore 
and we recover the conventional laser-matter formulation. The units of e 
depend on the function g{r). For instance, if ^'(r) = 2 ; (a linear function), e 
has units of inverse length (see e.g. 0 ). 

Often, it is desirable to solve the TDSE in the length gauge or the maximal 
coupling gauge. This is so because the numerical or analytical calculation 
can be expressed in an easy way and the computation of certain observables 
is more efficient [29]. Therefore, the main question is how we can perform 
the transformation of the Hamiltonian in the VG, Eq. ([^, to the LG. The 
gauge transformation should be boiling down in an unitary translation of the 
whole wavefunction [21|. We dehne this unitary transformation according to: 


Tv = 


( 4 ) 


where, Ty = Tv(r,f) and Tl = TL(r,t) are the wavefunctions in the VG 
and LG, respectively. Q is the unitary hermitian operator dehned according 
to the following rule Q = exp [fx(r, t)] [211 [32l |33|, with x(r, t) = A(r', t) ■ 

dr'. The latter expression is a contour integral which is independent of the 
path, because we can safely assume that the effect of the magnetic held is 
negligible, i.e. that the curl of the vector potential for the inhomogeneous 
held is zero, V x A = 0. Furthermore, by using Eqs. ([^ and (|^, we hnd the 
transformation for the Hamiltonian, Hy, from the VG to the LG: 


ot at 


( 5 ) 


Then, knowing that E(r,t) = —^A(r,f), i.e. the relationship between 


5 




E(r,t) and A(r,t), the last expression becomes: 

QHyQ^ + [ E{r',t)-dr' 

Jc 

As the TDSE is gauge invariant, we infer that the Hamiltonian in the VG, 
Hy, is transformed to the LG, Hi^, via: 

i/L = QHyQ^+ [ E(r',f)-dr'. (7) 

Jc 

It can be demonstrated that the hrst term on the right hand side of Eq. Q, 
yields QHyQ^ = + H(r)- Then, the Hamiltonian in the LG takes the 

form: 

Hl = ^ + Ho(r) + Gnt(r,t), (8) 

here, p is the kinetic momentum operator, and I4nt(r,t) = /cE(r',f) ■ dr', 
is a contour integral. In terms of Glassical Mechanics, we can interpret this 
last term as the work done in the electric held E(r',f) to move the electron 
from an arbitrary place to the position r. In the particular case when the 
vector potential has the functional form given by Eq. (|^ and the function 
g{r), is set to ^'(r) = z, the Hamiltonian operator in the LG becomes: 

Hh = ^ + To(r) + ^(1 +-2;)Eh(t), (9) 

where, Eh(t) = —denotes the spatial homogeneous part of the laser 
electric held. Gommonly this held, Eh(t), is called the conventional or spatial 
homogeneous held. 

In the next section, we shall compare the numerical accuracy of the VG 
and LG predictions for the high-order harmonic generation (HHG) driven by 
plasmonic helds. Our numerical models are based on Eqs. (|^, for the VG, 
and 0 . for the LG, respectively. 

3. Numerical algorithms 

The methods utilized to numerically integrate the TDSE are classihed 
by considering how the time evolution of the EWE is computed. When the 
EWE at a later time is obtained from the one at the current time, we have 


Tt = i 




( 6 ) 


6 





the so-called explicit methods. On the other hand, implicit schemes hnd 
the EWF by solving an equation involving both the actual EWF and one 
at later time. We choose the Spectral-Split Operator (SO) method joint 
with the Crank-Nicolson (ON) scheme, which are both explicit methods, to 
numerically integrate the TDSE of our interest. The SO uses a spectral 
technique to evaluate the derivative operator in the Fourier domain mm, 
and, on the other hand, the ON is based on the hnite element difference 
discretization technique |3l] to implement the second derivative present in 
the Hamiltonian, which dehnes the kinetic operator term. 

In order to test the accuracy of both the VG and the LG in the HHG driven 
by plasmonic helds, we have implemented the TDSE via the SO and GN 
techniques within a one spatial dimension model (ID). 

A general solution of the TDSE is done by employing a unitary U{to + At, to) 
evolution operator, where to is the initial time, i.e. the initial EWF \ko(to) is 
known and we evolve the system to an unknown state T {to + At) at a given 
time to + At [55] : 


T(to + Af) — f/(to + At, to)ft'o(^o)- (10) 


For simplicity, in Eq. (10), we have dropped out the spatial (r) dependence 


on the EWF. In the laser-matter community, the U{to + At, to) is commonly 
known as a propagator and it has the following explicit form, U{to + At, to) = 


exp 




/to 


In reference [5T], Feit et ah have introduced the SO method to numerically 
solve the TDSE in two spatial dimensions (2D) by using Eq. (10). This 


method consists in to split the time evolution operator U{to + At, to) 
three parts [5T] : 


^{to + At) = 


Here, the Hamiltonian, H{to + ^), is divided in H{to + ^) = + I4ff(0 + 

^), with Vef[{t) = Vb(i') + /c E(r', t)-dr' the effective potential in the LG. The 
main advantage of Eq. 10 is that we can evaluate the kinetic operator term, 
to), acting on the initial state, in the momentum space. This 
means that we need to compute a Forward Fourier Transform (EFT) [26] of 
'Lo(r, to) and then multiply it by a phase factor which evaluates the action of 
the kinetic operator, instead of a complicated derivate operator. Then, over 
this momentum space EWF, an Inverse Fourier Transform (IFT) is applied 
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in order to return to the coordinate space [SI]. This procedure is performed 
because the momentum operator in the conjugate (momentum) space is just 
a number and not a derivative one. 

For the conventional or homogeneous helds case, the vector potential does 
not depend on the spatial coordinate, i.e. A(r,t) = A(t), allowing us to 
evaluate the kinetic operator in the VG as: 

^{to + At) = e-*5{P+A(*o+At/2))2Ai/2g-iyoAtg-ii(p+A(to+At/2))^At/2^^j'^^^_ j']^2) 


Clearly, this is not the case for the spatial nonhomogeneous helds. The 
dependence of the vector potential on the position, as stated in Eq. ([^, does 
not allow us to apply Eq. (12). This is so because in the momentum space 


the position operator becomes a derivative, which complicates substantially 
the SO method. Therefore, we conclude that the SO method can not be 
easily employed to numerically integrate the TDSE in the VG. However, by 
using a hnite element grid discretization, we will show that the CN method 
can be used in both gauges, VG and LG. The CN method is based on the 


solution of Eq. (10) by the Caley formula and the evaluation of the kinetic 
operator in the position space using a hnite element method [M]. In ID the 
numerical algorithm can be written as: 


At 

1 + i—^H(to + At/2) 


T(to + At) = 


At 

i—H(tQ + At/2) 


^o(^o)- (13) 


The unknown EWE, \k(fo + At), is then computed by solving a tridiagonal 
system of equations. 


4. System description and results 

In Attosecond Science, high-order harmonic generation (HHG) is one 
of the most important phenomena. For instance, it is possible to synthe¬ 
size attosecond pulses or to obtain structural information about the atomic 
of molecular systems [1] from the HHG spectra. Therefore, we chose here 
this observable driven by conventional (homogeneous) and non-homogeneous 
helds to compare the accuracy of our VG and LG implementations. 

For simplicity, we restrict ourselves to a one dimensional (ID) model, 
although it is known this approach is able to accurately reproduce the main 
features of the HHG spectra of real atoms [36] • The potential well, Vo{x), 








which defines our atomic system, is a soft-core or quasi Coulomb potential: 


V,{x) 


Z 

+ a 


(14) 


where Z is the atomic charge and a a parameter which allows us to tune the 
ionization potential of the atom of interest. In this paper, we set Z = 1 and 
a = 0.488 a.u., such as the ionization potential is Ip = 0.9 a.u. (24.6 eV), i.e 
the value for the single active electron (SAE) model of the He atom [27]. Our 
ground state was computed via imaginary time propagation for a different 
set of spatial grid steps 5x. To assure a “good time” step, convergence, 
we have used the criterion: 5t < 5x^/2 (for more details see e.g. [34]b 

In order to compute the HHG spectra, we firstly calculate the dipole 
acceleration expectation value, as a function of time: 


+ (15) 

where the EWF T(x,t) is obtained via the SO and CN methods already 
described in the previous Section. The spectral intensity, /hhg(i^) = |od(i^)p! 
of the harmonic emission is then computed by Fourier transforming the dipole 
acceleration by using: 

/ + 00 

. (16) 

■OO 


The numerical computation of the HHG spectra will be performed by using a 
set of position steps = {0.05, 0.1, 0.15,...} a.u. Consequently, and in order 
to estimate the numerical convergence of the HHG spectra as a function of 6x, 
we use the spectral intensity difference between the smallest step, i.e. 5xq = 
0.05 a.u., and the rest of the set, = |7hhg,5xo(i^) - 7HHG,5xj(t<;)|, 

with j = (1,2,...}. Furthermore, for each of the 6x, the computing time is 
also retrieved for both the VG and LG. 


4-1. HHG driven by conventional fields 

Firstly, we present computations of the harmonic spectra intensity, /hhg(i^), 
driven by a conventional homogeneous field. This case is the limit e —)■ 0. 
Both numerical methods above described, i.e. the SO and CN, have been 
used to compute the emitted harmonic spectra intensity both in the VG and 
LG. We shall show below that both gauges give the same results. 
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Figure 1: (color online) Computed high-order harmonic intensity spectra (in arbitrary 
units) driven by a spatial homogeneous (conventional) laser field under the LG (red solid 
line) and the VG (blue dashed line). Panel (a) the HHG spectra are obtained by using 
the SO method, panel (b) the same as (a) but using the GN method. The green vertical 
dashed line depicts the classical harmonic cut-off law, i.e., ric = {Ip+ 3.17 Up)/ ujq [35]. The 
laser pulse parameters for these simulations are: intensity Iq = 2 x 10^'* W/cm^, carrier 
frequency ujq = 0.057 a.u. (corresponding to a wavelength of A = 800 nm), and GEP, 
7>CEP = 0 rad. The pulse envelope is a sin^ function with four total cycles. We chose a 
grid step of Sx = 0.05 a.u. for both gauges. 


The TDSE calculations are performed in a grid with a step 6x = 0.05 
a.u., and a spatial grid length of = 3500 a.u. The real-time evolution is 
done with a time-step of 6t = 0.00125 a.u. Fig. [TJshows the spectral intensity 
of the harmonic emission when a laser pulse interacts with our ID helium 
model. In Fig. 1(a), the comparison of the HHG spectra between the LG and 
VG is depicted by using the SO method. The same comparison is shown in 
Fig. 1(b), but here the GN method is used for the numerical integration of 
the TDSF. Both methods show a perfect agreement when the LG and VG 
are used to compute the spectral harmonic intensity. This confirms that our 
numerical methods are able to describe the HHG process for any of the grid 
steps used in our simulations. 

As a next test, we have integrated the TDSE for a set of grid steps, 6x = 
{0.05, 0.15, ...,0.45} a.u., and computed the emitted harmonics. Figure]^ 
shows the results of the harmonic intensity, /hhg(i^)! as a function of the grid 
step computed by the SO Figs. |^a)-(b) and the GN Figs. |^c)-(d) methods, 
respectively. For both the VG and LG, the numerical spectra by using the 
SO method shows a perfect agreement for all the set of grid steps, Sx used in 
our simulations. In contrast, the situation is different when the GN method is 
employed to compute the harmonic spectra. For instance, the Figs. c)-(d) 
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Figure 2: (color online) Computed high-order harmonic intensity spectra driven by a 
spatial homogeneous (conventional) field (in arbitrary units) as a function of the grid 
step, Sx, under the LG and the VG by the SO method, panels (a)-(b) and the CN method 
panels (c)-(d). The green vertical dashed line depicts the classical harmonic cut-off law. 
The laser pulse parameters for these simulations are the same as in Fig. 1, i.e. intensity 
/q = 2 X 10^"^ W/cm^, carrier frequency Wq = 0.057 a.u. (corresponding to a wavelength 
of A = 800 nm), and CEP, i^cep = 0 rad. The pulse envelope is a sin^ function with four 
total cycles. 
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show that the emitted harmonic spectrum depends on the grid step when the 
LG and VG are employed to computed the HHG. In addition, the computed 
HHG spectra slightly differ in the whole harmonic-order range whether the 
LG or the VG is used in the calculation and for the larger grid steps, i.e., 6x > 
0.25 a.u. Additional structures can be observed in the low-order harmonics 
for the case of VG (see Figs. id)). although the general shape, including 
the harmonic cutoff, is in excellent agreement with the rest of the schemes. 
Gonsidering the numerical error that the finite element method has for the 
second derivative as a function of the grid spacing 6x, it is reasonable to 
attribute poor convergence when the GN method is used with larger grid 
steps 6x. Furthermore, in view of the fact that the VG has an extra spatial 
derivative of first order within the Hamiltonian, p ■ A, we would expect that 
the numerical accuracy decreases when the spatial step, 6x, increases. This is 
the reason behind the noticeable difference between the LG and the VG when 
larger grid steps are employed in the calculations of the HHG. Our numerical 
results show, however, that this difference between LG and VG disappears 
for the smallest spatial grid steps, i.e., 6x < 0.2 a.u. These outcomes suggest 
that the best method to compute the HHG spectrum is the SO. On the other 
hand, in cases where the SO method is challenging, the GN method can be 
used if the grid step is small enough, e.g. 5a: < 0.1 a.u. We should note that 
the grid step will depend on the particular problem, i.e. laser parameters, 
etc., although we can expect a general trend. For this reason, we suggest to 
perform a convergence analysis if the GN method is employed and to chose 
the adequate parameters for the required accuracy. 

In the next, we shall perform the computation of the HHG spectra driven 
by a spatial inhomogeneous field. For the reasons explained in Section 3, we 
shall only use the GN method and compute the harmonic emission both in 
the LG and VG. 

4-2. HHG driven by spatial inhomogeneous fields 

As was mentioned at the outset, when a laser field is focused on a metal¬ 
lic nanostructure, a hot spot of higher intensity, high enough to exceed the 
threshold for HHG in atoms, is created due to the coupling between the in¬ 
coming field and the surface plasmon polaritons (SPPs) [2]. The main prop¬ 
erty of the effective laser electric field is that it presents a spatial variation in 
the same scale as the one of the dynamics of the active electron. Therefore, 
the interaction between this plamonic field and the atomic electron, which 
governs the HHG process, will change substantially. As the electric field is 
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Figure 3: (color online) HHG spectra driven by a spatial inhomogeneous field computed 
by using the LG (blue line) and VG (red dashed line). The parameters for the laser pulse 
are the same that those used in Fig. the inhomogeneous parameter is e = 0.0175 a.u. 
(see the text for more details) and the grid step is 6x = 0.05 a.u. 


not anymore spatially homogeneous, the electron will experience different 
electric field strengths along its trajectory. The question that emerges is 
which gauge can give us a numerical advantage when the TDSE is solved for 
the computation of the HHG spectra driven by spatial inhomogeneous fields. 
Before to address this question, we firstly demonstrate that both the LG and 
VG are equivalent in the description of the HHG driven by nonhomogeneous 
fields, as was demonstrated by the conventional case (see Section 4.1). 

We have numerically integrated the TDSE in ID for the same atomic sys¬ 
tem used in the previous section (Section 4.1), but now the effective electric 
field is spatially inhomogeneous. Fig. shows the comparison between the 
calculated HHG spectra driven by an inhomogeneous field for both the LG 
and VG. The inhomogeneous parameter value is set at e = 0.0175 a.u., which 
correspond to an inhomogeneous region of about 60 a.u. (3 nm) (see [5] for 
more details). Perfect agreement between the predictions of both the LG 
and VG are found. Therefore, these results suggest that our derivations are 
appropriate for spatial nonhomogeneous fields as well. As a consequence, 
this invariance allows us to check which gauge can be more convenient to 
compute the HHG driven by spatial nonhomogeneous fields. We will address 
this point by considering the convergence of both the LG and VG. In other 
words, which of the two gauges presents less numerical error against the grid 
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Figure 4: (color online) HHG spectra driven by plasmonic fields both in the LG and 
VG as a function of the grid step Sx are depicted in panel (a) and (b), respectively. We 
have used the CN method to numerically integrate the TDSE. The parameters for the 
laser pulse are the same that those used in Fig. and the inhomogeneous parameter is 
e = 0.0175 a.u. 


step, 6x, and which one is faster in the compntation of the HHG spectra. 
Fig.g shows the HHG spectra as a fnnction of the grid step for both the 
LG Fig. 4(a) and the VG Fig. 4(b) compnted by nsing the GN method. The 
HHG spectra for the LG show a convergence for the smallest grid step, i.e., 
for 6x = 0.05 a.n. We shonld note, however, that the highest freqnency 
of the HHG spectra change when the grid step is increases, which suggests 
that the computation of the HHG spectra driven by spatial inhomogeneous 
fields deserves special attention when “large” grid steps are used. A similar 
result is found when the VG is employed although it is possible to observe 
convergence for larger values of 6x. A suitable way to conhrm the HHG 
cutoff and corroborate the convergence of the numerical schemes, is to rely on 
classical simulations. It is known that the limits on the HHG spectra can be 
obtained via classical simulations, e.g. by computing the maximum electron 
kinetic energy upon recombination |3S]. For spatial nonhomogeneous helds, 
it was demonstrated a perfect agreement between the classical predictions 
and the TDSE simulations (see e.g. 0 ) and, as a consequence, we could 
benchmark our VG and LG approaches by solving the classical equations of 
motion for an electron moving in an oscillating and spatial dependent electric 
field (for more details see [25]). 

In addition, despite of the fact that for larger grid steps the LG shows 
a large deviation for the highest frequency compared to the VG results, we 


evaluate the relative error defined by 


^HHG,c5xq 


for each gauge as a function 
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Figure 5: (color online) (a) Convergence relative-error as a function of the grid step by 
using the LG (blue line with squares) and the VG (red line with circles), (b) computing 
real-time as a function of the grid step for both the LG (blue squares) and VG (red circles). 
The simulation parameters are the same that those used in Fig. 


of the grid step 5x. The results are depicted in Fig. |^a). This panel shows 
that a large difference appears whether LG or VG is used to compute the 
HHG spectra by the GN method. For values of 5x larger than 0.2 a.u., the 
relative error between the LG and the VG has a difference of about two 
orders of magnitude, which suggests that the LG would be more appropriate 
than the VG to compute the HHG spectra. Finally, in the panel (b), we show 
the computational time for each gauge. As can be observed the computing 
times for both the LG and the VG are similar. From this consideration we 
can conclude that the LG could be the most appropriated gauge in order to 
compute the HHG, given the fact it allows us to use larger grid steps. 

5. Conclusions 

We have reviewed the gauge invariance problem, both analytical and nu¬ 
merically, for the calculation of the HHG phenomenon driven by spatial ho¬ 
mogeneous and inhomogeneous (plasmonic enhanced) electric helds. To this 
purpose we have solved the TDSE in reduced dimensions by implementing 
the Spectral-Split Operator and the Grank-Nicolson algorithms. It was found 
that both the LG and VG are equivalent in the description of the harmonic 
emission processes for each of the two studied cases: the spatial homogeneous 
and the spatial inhomogeneous helds. For the spatial inhomogeneous held 
case, and due to the dependence of the vector potential on the position, we 
found that the SO method was difficult to implement in the numerical solu- 
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tion of the TDSE. In contrast, the CN method has shown advantages because 
it is based on a hnite element discretization. Our numerical results based on 
the CN method suggested that the calculation of the harmonic spectra de¬ 
pends strongly on the grid step chosen to perform the numerical integration. 
Both gauges are equivalent, but according to the numerical convergence of 
the HHG spectra, the LG apparently appears to be more accurate than the 
VG for the lowest harmonics. This is so because the lowest harmonics change 
by several orders of magnitude when the grid step increases. Furthermore, 
it was shown that particular attention in the choice of the spatial grid step 
should to be taken when spatial inhomogeneous helds are employed. This is 
so, because the limits of the harmonic radiation appear to be very sensitive 
to this parameter. Finally, it was found that the computational time was 
similar for both the LG or VG, if they were used for the computation of the 
HHG spectrum in the moderate and high laser intensity regimes. 
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